rm(list = ls())

library(microsynth)
library(xtable)
library(tidyverse)

rank_fun <- function(x){
  su = sort(unique(x))
  for (i in 1:length(su)) x[x==su[i]] = i
  return(x)
}

## Load Data 

dat <- readRDS('data/nrw_rlp_county_1990_2021_clean.rds') %>%
  filter(substr(county_id, 1, 2) == '07') %>%
  mutate(period = rank_fun(year)) %>%
  mutate(greens_zweit = greens_zweit*100) ## outcome in p.p.

## Drop weakly treated areas 

dat <- dat %>%
  filter(treat_categ %in% c('none', 'sehr schwer'))

## Define binary treatment 

dat <- dat %>%
  mutate(treat_very_heavy = ifelse(county_id == '07131' & year >= 2021, 1, 0))

## Synthethic control 

sc_rlp <- microsynth(dat, 
                     idvar= 'county_id', 
                     timevar = 'period', 
                     intvar = 'treat_very_heavy', 
                     result.var = c('greens_zweit'), 
                     start.pre = 1, 
                     end.pre = 8, 
                     perm = 50,
                     check.feas = TRUE, 
                     use.backup = TRUE,
                     jack = F, 
                     test = "twosided")


plot_dat <- sc_rlp$Plot.Stats 

## Weighted mean in treated/control groups

synth_treated <- plot_dat$Treatment %>%
  t() %>%
  as_tibble() %>%
  mutate(group = 'Ahrweiler county',
         period = unique(dat$year))

synth_control <- plot_dat$Control %>%
  t() %>%
  as_tibble() %>%
  mutate(group = 'Synthetic control',
         period = unique(dat$year))

## Combine 

plot_no_se <- bind_rows(synth_treated, synth_control)

## FIGURE 5

p1 <- ggplot(plot_no_se, aes(x = as.numeric(period), y = greens_zweit, 
                             group = group, col = group,
                             fill = group)) + 
  geom_line() + 
  geom_point(shape = 21) + 
  theme_bw() + 
  labs(x = '',
       y = 'Green party vote share (%)') + 
  theme(legend.title = element_blank(),
        legend.position = 'bottom') + 
  scale_color_manual(values = c('black', 'gray')) + 
  scale_y_continuous(limits = c(0, 14), breaks = seq(0, 14, 2)) + 
  scale_fill_manual(values = c('black', 'gray')) + 
  geom_vline(xintercept = 2020, linetype = 'dotted') + 
  scale_x_continuous(breaks = unique(as.numeric(dat$year)))

p1

## TABLE A 15 

weights_sc <- sc_rlp$w$Weights[,1] %>%
  as_data_frame() %>%
  mutate(county_id = row.names(sc_rlp$w$Weights)) %>%
  rename(weight = value) %>%
  mutate(weight = round(weight, 4)) %>%
  filter(county_id != '07131') %>%
  arrange(-weight) %>%
  left_join(readRDS('data/counties_list.rds')) %>%
  dplyr::select(county, county_id, weight) 

weights_sc %>%
  xtable::xtable()




